Magnetic phase diagram of a spin-1 condensate in two dimensions with dipole 

interaction 
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Several new features arise in the ground-state phase diagram of a spin-1 condensate trapped in 
an optical trap when the magnetic dipole interaction between the atoms is taken into account along 
with confinement and spin precession. The boundaries between the regions of ferromagnetic and 
polar phases move as the dipole strength is varied and the ferromagnetic phases can be modulated. 
The magnetization of the ferromagnetic phase perpendicular to the field becomes modulated as 
a helix winding around the magnetic field direction, with a wavelength inversely proportional to 
the dipole strength. This modulation should be observable for current experimental parameters in 
*^Rb. Hence the much-sought supersolid state, with broken continuous translation invariance in 
one direction and broken global U{1) invariance, occurs generically as a metastable state in this 
system as a result of dipole interaction. The ferromagnetic state parallel to the applied magnetic 
field becomes striped in a finite system at strong dipolar coupling. 

PACS numbers: 03.75.Hh, 03.75.Mn, 03.75.Lm 



I. INTRODUCTION 

Bose condensates of atoms with nonzero total spin 
F > 1 show various phases combining magnetic and su- 
perfluid order. When the magnetic symmetry is broken 
spontaneously, as can occur when the atoms are confined 
in a spin-independent optical trap, condensates are clas- 
sified as "polar" (for antiferromagnetic interactions) or 
"ferromagnetic" . Most theoretical studies of these spinor 
Bose condensates neglect the long-range interaction be- 
tween atomic magnetic moments, and this neglect is jus- 
tified for many experimental conditions. However, re- 
cent experimcnts^^'^ investigating ordering in a nearly 
two-dimensional condensate have shown complex mag- 
netic behavior in the ferromagnetically interacting F ^ 1 
spinor Bose gas of *^Rb. 

The most surprising feature of these experiments, 
which image the spin distribution in real space, is a long- 
lived phase that appears to have the broken global t/(l) 
invariance of a superfluid along with possible breaking 
of the continuous translational symmetry in one or 
two directions, i.e., with stripe-like or checkerboard- like 
order. A possible supersolid phase has recently also been 
suggested in the superfluid of ^He,^ Many theoretical 
papers have been written about the properties of ''He 
and whether a supersolid phase can exist in the absence 
of disorder. Only recently have theoretical studies been 
done to explain the observed supersolid-like behavior in 
a *^Rb spinor condensate.^ The earlier studies of *^Rb 
concentrated on magnetic properties arising from the 
weak spin-dependent local interaction and the quadratic 
Zeeman shift. More recent experiments^"^ indicates that 
the long range dipole interaction also plays an important 
role in the formation of the magnetic phases in spatially 
large systems and with this addition a supersolid state 
might be possible. 



Most previous studies of this system concern dy- 
namical properties: the leading instability when the 
Hamiltonian is changed to favor ferromagnetic order 
can be stripe-like or checkerboard-like depending on 
parametersi^"— In this paper, our goal is to determine 
the static ground-state phase diagram. We start from 
the phases that are well established at low tempera- 
tures^""— for a spin-1 gas with no dipole interaction 
and quadratic Zeeman effect. (Low temperatures mean 
below the superfluid and magnetic transitions, where all 
the studies in this paper will take place.) We then add 
the dipole interaction to see how it changes the phases 
as well as the location of the boundary between them. 
We do this in a quasi-two-dimensional geometry as in 
the experiments.^"'^ We investigate both an infinite and 
a finite square planar geometry. After observing the 
formation of two kinds of stripe order in a Monte Carlo 
simulation, we developed an analytical approach to 
explain the results, based upon smallness of the dipolar 
coupling at short distances. That analytical approach is 
presented first in order to prepare the groundwork for 
the simulation results. 

We show that all boundaries in the phase diagram, ex- 
cept between the two polar phases, are moved when the 
dipole interaction is added, some in a non-intuitive way. 
The magnetic dipole interaction prefers a ferromagnetic 
state, but the confinement makes a ferromagnetic state 
out of the plane energetically unfavorable. Moreover, the 
spin precession make the in-plane perpendicular ferro- 
magnetic state unfavorable, since the spin rotates out of 
the plane. Both ferromagnetic phases can get modulated 
in one direction. The phase parallel to the external fields 
needs a strong dipole interaction or a system much wider 
than its length to become modulated. This modulation 
appears as fully magnetized stripes with sharp domain 
walls between them. The phase perpendicular to the ex- 
ternal fields gets modulated, from the very lowest dipole 
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strengths, into a helical configuration around the field. 
The wavelength of the helix is inversely proportional to 
the dipole strength. For *^Rb the wavelength is ~ 80/im 
and should be observable in experiments. 

The outline of this paper is as follows. In the following 
section, we review the basic physics of spinor conden- 
sates without the dipole interaction. In Section III we 
introduce the dipole interaction and put it into a form 
that is convenient for numerical simulations. Section IV 
presents analytical results in the limit of weak dipole in- 
teraction, and Section V contains the results of our Monte 
Carlo simulations of the problem. The final section sum- 
marizes the relationship between our results and those of 
other theoretical papers and suggests how future experi- 
ments could be designed to observe clearly the metastable 
supersolid phase found in our simulations. 



II. REVIEW OF SPINOR CONDENSATE 
WITHOUT MAGNETIC DIPOLE INTERACTION 

A Bose-Einstcin condensate of spin F — 1 atoms is de- 
scribed by a three-component complex order parameter 



*(x) = v/?^3d(x)-0(x) = Vn3D(x) 



^+i(x) 
V'o(x) 
V^-i(x) 



(1) 



where the spinor ■(/'(x) is normalized as tp^ip = 1 and the 
subscripts label the spin eigenvalue with respect to an 
arbitrarily chosen quantization direction. In the absence 
of external fields and neglecting the dipole interaction, 
the Hamiltonian governing the condensate 



Ho 



Co 



'3D 



, (2) 



where m is the atomic mass, M(x) = '(/'^(x)F'0(x) is the 
dimensionless magnetization (|M| < 1), and {f } are the 
three generators of SU{2) in the spin-1 representation 




(3) 



The first term in the Hamiltonian is the kinetic en- 
ergy for bosons with mass m. The next two terms 
are the spin-independent and spin-dependent contact in- 
teractions, respectively. The coefficients are given by 
Co = (47rfi^/3m)(2a2 -I- ag) and C2 = {Anh? /3m){a2 — uq), 
with {oq, 02} the s-wave scattering lengths in the channel 
with total angular momentum {0, 2}. 

When C2 < ("ferromagnetic") it is energetically fa- 
vorable for this system to magnetize, M ^ 0, while 
C2 > favors a "polar" state with M = 0. The scattering 



lengths for ^'^Rb are Uq = 101. Sas and 02 = lOOAasr^ 
where as is the Bohr radius, so C2 is negative and its con- 
densate will be ferromagnetic in the absence of external 
fields (still neglecting the dipole interaction). However, 
the condensate of ^"^Na will be in a polar state^S 

The external fields normally applied to a spinor con- 
densate consist of an optical trap and a uniform magnetic 
field described by the following addition to the Hamilto- 
nian 



Hef^ Id 



(4) 



The trapping potential U (x) confines the condensate spa- 
tially; for our purposes, its main effect will be to produce 
a quasi-two-dimensional geometry. The quadratic Zee- 
man shift q can be tuned independently of B with mi- 
crowave radiation, q — + q^^^ }^ We take the two 
sources as coaxial along z, so we can use Eq. (U). This is 
also the axis we quantize the spinor along. The magnetic 
field also creates a linear Zeeman term B • J(Px nsjj^, 
that favors an uniformly magnetized condensate. How- 
ever, experiments on ^^Rb have not observed any ten- 
dency toward such relaxation over the accessible time 
scales of several seconds making the longitudinal com- 
ponent of magnetization conserved. (This assumption 
does not apply in condensates of higher spin, such as 
chromiumJ^) Normally, this component is chosen to van- 
ish initially and can hence be ignored for the purpose of 
energetics. However, the magnetic field also causes Lar- 
mor precession of the magnetization perpendicular to it. 
This is an important effect that needs to be taken into 
account as it modifies the nature of the magnetic inter- 
action on time scales longer than the precession time. 
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FIG. 1: (Color online) Examples of possible spin configura- 
tions in the plane. The external fields are along the horizontal 
axis, Mzix, z) is plotted on the horizontal axis, Mx{x, 2) is on 
the vertical for every plaquette and My{x, z) is not shown, (a) 
Uniform fully magnetized F\\ , (b) striped fully magnetized , 
(c) uniform partly magnetized F±/P\\ state, (d) helical fully 
magnetized Fx • 
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The spin state of the condensate is governed by the 
parameters C2 and q, as in Fig. There are two dif- 
ferent kinds of polar states (c2 > 0), one that minimizes 
{{F^)^) = and one that maximizes {{F^)"^) = 1 the 
impact of the quadratic Zeeman term. Respectively, 





(5) 



Consequently, does the phase , with order-parameter 
manifold C/(l), appears at g > 0, while the phase P± 
appears when q < 0. Note that the range of 9 is only 
[0, tt), or alternatively that the order-parameter manifold 
for this phase is U{1) x [/(1)/Z2.^5- When C2 < and g < 
0, both energies are minimized by ferromagnetic states 




(6) 



giving a manifold U{\) x Z2 (recall that we exclude the 
linear Zccman energy from energetic considerations), see 
Fig. [1] In the final quadrant of the phase diagram, how- 
ever, no ferromagnetic state minimizes the quadratic Zee- 
man energy. The smallest impact of a ferromagnetic state 
on the quadratic Zeeman term is ((P^)^) = 1/2 for 



V2 



(7) 



Consequently, for q < qc = 2|c2|n3D the state will be 
a linear combination of i/'y' and 'ipA_ with magnetization 

M^ + iMy = - {q/qcY and manifold U{1) x [/(I), 
see Fig. [TJ Above q^ the state wiU be the pure polar 
state i^i . 

Typical experimental values for ^^Rb^""^ include a peak 
density of no = 2.5 x lO^^'^cm"'^, giving the interaction 
strengths cotlq ~ 1.9 kHz and C2no = —9 Hz, while q^ « 
1.6Hz and q^^^ can be tuned from roughly —50 Hz to 
50 Hz and is normally taken coaxial to q^ M 




FIG. 2: (Color online) Ground state phase diagram of a spin- 
1 condensate without dipolar interaction; from Mukerjee et 



profile and a Gaussian, 

Piiy) - ^Q{K/2-y)Q{K/2+y), p^iy) = —^[^e~'^, 

A ay \ TT 

(9) 

where 8(x) is the Heaviside step function. We introduce 
a common notation for the condensate thickness T and 
a 3-dimcnsional density n^]:){x,z) without y-dependence 

i = {p)= Jdypiyf, nM^,z)^!^^£p^ (10) 

for the boxcar profile and for the gaussian profile, to be 
able to treat both profiles simultaneously in section ITVl 
In most of our analysis these densities are also indepen- 
dent of {x,z), except where we use a nonzero trapping 
potential U{x, z) in the plane. 



B. Precession 



A. Confinement 

The optical trap in the experiment makes the gas ef- 
fectively two dimensional, with a Thomas-Fermi radius 
ttf ~ 1.5pm along the direction of tightest confine- 
ment i^Ti^ Since this is smaller than the spin healing length 
^ = ■y//i^/(2m|c2|n3D) « 2.5pm, we take the confinement 
to be along the y direction and treat the gas as frozen 
along this direction; that is, we take 

*(x) = ^/n2D{x,z)p{y)ip{x,z), (8) 

where we assume Jdy p{y) = 1. In the following we will 
consider one of two profiles p{y) as convenient, a boxcar 



Atoms with magnetic moment fi± = gppsMj^ per- 
pendicular to the field precess at frequency I7I-B0 = 
IgplfJ-sBQ around the fields. As usual, pb is the Bohr 
magneton and gp is Lande's g-factor. For ®^Rb, gp = 
— 1/2 and a field of Bq = 150 mG produces a Larmor 
precession at 110 kHz, a scale orders of magnitude larger 
than the contact interactions or the quadratic Zeeman 
energy. 

The Hamiltonian considered so far is invariant under 
the spin rotation 

Mx,z)~^Uki{t)M^,z), C/(i) = e-»^S"B.Ft 

and is hence unaffected by the rapid Larmor precession. 
Therefore, adding precession does not affect the phase 



4 



diagram in the problem with only local interactions 
However when we include the dipole interaction in the 
next section, both confinement and spin precession be- 
come important. 

III. MAGNETIC DIPOLE INTERACTION 

The interactions considered thus far for a spin-1 con- 
densate are all local. However, the moments /x will inter- 
act through the long-ranged dipole interaction. This is 
weak for ®^Rb relative to most other energies in the sys- 
tem, but since it is long ranged it will have an important 
impact on the magnetic phases. The initial studies of the 
spin-1 condensate ignored this term,^'^- but some recent 
works have included it along with the effects of quasi- 
two-dimensional confinement and rapid Larmor preces- 
sion. ^'^ Among other results, it was shown that dipolar 
interaction renders the Larmor precession unstable,^ and 
we return to this point in the concluding section. Un- 
til then we follow previous authors and assume that this 
instability has significant effects only at late times, and 
so neglect it. Cherng and Demler examined the instabil- 
ity spectrum of a uniform ferromagnetic state within a 
mean field and collective mode analysis. We will use the 
same physical model but instead look at the ground state 
phase diagram and consider a wider range of parameters 
C2, q, and Cd (see Eq. ((T^ below) with analytical and 
Monte Carlo calculations. 

The total Hamiltonian we work with is 



IV. ANALYTICAL RESULTS 



H — Hq + Hef + Hd 



(12) 



where 



Hd^p = ^J d^xd^x'n3D{x)M,{x)n3D{^')Mj{x') 



f%<5(3)(x-x') 



(13) 



This is the same as the more usual expression with 
{Sij — 3'rifj)/r^, but split it into a part that is positive- 
(semi-) definite and a part that simply shifts the parame- 
ter C2 — C2 — 47rC(i/3 (see the beginning of Appendix Elfor 
a fuller discussion of the magnetic dipole term). Indeed, 
with two integrations by parts the first term becomes the 
Coulomb interaction for a charge density V • (ria/^M). 
We will typically mean just this term when referring to 
"the dipole interaction," since it is the difficult part. The 
strength of the dipole term is given by Cd = ^o5FMB/47r, 
where /io is the vacuum permeability, giving a value of 
Cdno = 0.8 Hz for *^Rb. The effect of confinement is 
less trivial for this term then for the others, and trans- 
forming to a rotating frame is also nontrivial since the 
interaction couples spin directions to spatial directions. 
See Appendix |^ for a full treatment of these effects. In 
the following section, we discuss how the dipole interac- 
tion is expected to modify the phase diagram when it is 
sufficiently weak that the Hdip = ground states can be 
used as a starting point. 



Adding the dipole interaction Eq. (|13l) will change the 
phase diagram Fig. The term that looks like the spin 
dependent interaction will just move the whole phase- 
diagram up along C2 with The energy from the 
Coulomb part of the dipole interaction is always positive, 
hence this parts prefers a polar state with zero magneti- 
zation M ~ 0. Consequently, regions of Fig.[2]with polar 
states above C2 = will not change if we add the 
dipole-dipole coupling. However, the rest of the phase 
diagram may be affected and the phase boundaries will 
depend on Cd, as we now discuss in some detail. 



A. Weak dipole interaction 

Adding a weak dipole term (weak compared to the ki- 
netic energy term) will only change the phase diagram 
slightly. We start out by ignoring any new phases and 
investigate how a weak dipole interaction will move the 
boundaries between the existing phases. The three mag- 
netic terms in the Hamiltonian are the spin-dependent 
contact interaction, the quadratic Zeeman and the dipole 
term. By comparing the energy contributions from these 
three for simple Ansdtze we can locate the boundaries 
between different minima, in a system with L the extent 
along z and W the extent along x. 

The polar phases are, of course, the simplest (see 
Eq. ©) 



£;^il = 0, E^^ = qn3DLWT. 



(14) 



Consider next the phase i^y , which appeared at g, C2 < 
in the system without dipolar energy. The effective 
charge density for such a state describes two quasi 1- 
dimensional lines of charge located at ±L/2, of length 
W. The self-energy of such lines of charge is given by 
2cd{n3DMfWT^ In W/T, see AppendixO to leading or- 
der. The other two terms are easily kept exact. Keep- 
ing terms of order and Aln^ where A = L, W, (see 
Eq. ©) 

= ^nlj^LWT + qnsDLWT, +2cdnljyWT^ \nW/T 

(15) 

with C2 = C2 - incd/S. 

The transition in the left half-plane betwen the states 
and P±, see Fig.|31 will hence be moved up from C2 = 

for a system without dipole interaction to 



/TT 

C2c = 4Crf y- - EL 



(16) 



where = ^'^^^ will vanish in the large-system limit. 

The region of the phase diagram with g > and C2 < 0, 
is the most interesting, due to the rapid precession of the 
perpendicular magnetization about the magnetic field, 
and the high dipolar energy cost of spins pointing out 
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FIG. 3: (Color online) Ground state phase diagram for a spin- 
1 condensate with dipole interaction and external fields, that 
introduces a quadratic Zeeman term and rapid spin preces- 
sion. Both polar and ferromagnetic phases appear, perpen- 
dicular as well as parallel to the field. 



of the plane. Consequently, the region of F^/ P\\ in the 
phase diagram will shrink and the regions of P\\ and F\\ 
grow, with the latter extending to positive values of q. 
For a uniform condensate with spins out of the plane, 
the Coulomb energy is equivalent to that of a parallel- 
plate capacitor, giving an energy 27rcd(ri3M)^T^LM^/T 
to leading order, i.e., neglecting fringing fields, see Ap- 
pendix [Cj 

Because of the precession, the spins will effectively 
average the out-of-plane and in-plane interaction ener- 
gies with equal weights. Consequently, the dipole energy 
for magnetization perpendicular to the external fields is 
Cd{n^Mf{-KW + T\nL/T)LT. To find the energy for 
the F±/P^\ state, we first have to find M, since this 
state is not completly magnetized. Consider a spinor 
V'^ = (a, 6, a) with a = y/{l - b^)/2 (1/^2 < 6 < 1), 
which represents a superposition of and (see 



Eqs. ([5]) and (O). Its magnetization is 
Putting it all together. 



2bVl~¥. 



eF^/P\\ ^ Al?{l-h^){^W+Cd{iiW+T\nWlT))ns,DLT 

+ qn3D{l-b^)LWT. (17) 
The energy for this state is minimized at 

.o 1 



1 



(18) 



As the notation suggests, the transition between the 
phases P\\ and F^/P\\ occurs ai q = qc, where iJ^J-Z^ii — 
E^w = and M = 



qc = 2|c2|n3D - AcdTizD y-^ 



ew 



(19) 



where ew = ^'w/t ^^^^ vanish in the large system limit. 
As can be seen in Eqs. ([T^ . ([TO]) and (1^ . the value of 
the magnetization and hence the order parameter for the 
P±/P\\ state decreases continuously and is zero at the 
phase transition to the Pn state. 



Mo = \{FjPn)\ = ^l-iq/q^y 



(20) 



This is exactly the same equation as for a system with- 
out dipole interaction, except that qc now is given by 
Eq. ([111). 

Plugging the form for b, Eq. ^TE\\ and p^ . back in also 
allow us to locate the transition between F^/P\\ and Fy, 
where E^^/^^^ = E'^n , which will occur at 



qc2 



\jqc (2|c2|n3D -HScdfiaD - 2eL)) ~ qc. (21) 



Finally, the transition between F|| and P\\ will take place 
when E^w = E^w = 0, at 



(22) 



The three transition lines (qc, 9c2 and q^) separating the 
three phases in the lower right quadrant meet at the point 

C2) = Acd (n^D + - 2ei^ , (^^ + 2ew - el 



.(23) 

To finish the phase diagram, we see that the transition 
line in Eq. ([25]), that separates Py and can be ex- 
tended to the region q,C2 > 0, with the substitution 
|c2| — > — C2 and that it will intersect with the transition 
line in Eq. ([T6)) at the point (g, C2) = (0, C2c). 



B. Magnetization textures 

The dipolar energy favors spatially modulated fer- 
romagnetic states, which screen the long-ranged in- 
teraction, over uniform states. Consider the state 
We can adapt a classic argument of Kittel con- 
cerning the formation of magnetic domains to the 
present quasi-two-dimensional geometryJi The bound- 
ary energy 2cdn\]^WT^ \x\W/T from before will become 
2cdn^jjWT^\nd/T if the uniform state breaks up into 
Ising-like domains of width d and length L that alter- 
nate between — 1 and = — 1, keeping the total 
magnetization Mq = 1 everywhere, see Fig.[TJ There will 
be a cost in kinetic energy at the domain walls, and the 
competition between these two effects sets the domain 
size. 

We can estimate an upper bound for the domain wall 
energy by assuming its width is the spin-healing length 
^5. The energy will scale with the area of the wall ~ LT, 
and the surface density will be aw ^ fi^n3D/2m^s- With 
the number of domains given by W/d, the energy is 



E 



aw- 



LWT 



2cdnloWT'^ \nd/T, 



(24) 
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which gives 



At leading order in the dipole strength, then, 



-L. 



(25) 



1 CdU^pT 
^ ri?/2m 



(26) 



The resulting domains have a width proportional to the 
length of the system, and are very large when the dipo- 
lar coupling is weak. In ^^Rb with the experimental pa- 
rameters given in section [II Al aw ~ 10^Hz/im~^ and 
(ill ~ 20L, which could be difhcult to achieve experimen- 
tally. 

For a rectangular sample (L > W) in the Fy state, 
with a constraint of zero total longitudinal magnetiza- 
tion (J dxn3£)(x)M^(x)) = 0), it can be more energeti- 
cally favorable to split up into two domains perpendic- 
ular to the field. The energy for this configuration is 
E = 2awWT + 3cdnljjWT^ \nW/T to leading order and 
if this is lower than the energy in Eq. (PT|) it will occur. 
However, this is only due to the constraint; a domain- free 
configuration has lower energy and a configuration with 
several domain walls perpendicular to the field will not 
be favorable for any values in the phase diagram. 

For the F± state, a different modulation will appear. 
In particular, since the state is XY-like (the rapid Lar- 
mor precession gives the same energy for all perpendic- 
ular spin directions), it can adopt a smoothly varying 
magnetization texture. The smoothest form will be a he- 
lix, with wave vector along the magnetic field, see Fig. [T] 
In other words, as shown in Fig. |4l the magnetization 
will adopt a configuration like Mx{z) — sai{kzz) and 
My{z) = sm{kzZ+ ^) at any instant of time. The ki- 
netic energy of such a state goes as fc^, while the dipole 
energy turns out to decrease as kz for small fc^; see Ap- 
pendix |B] for details. 




FIG. 4: (Color online) Transverse magnetization as a function 
of z, from numerical simulation. Orange: magnitude of total 
magnetization A/o, blue: transverse magnetization in plane 
Ma; and black: transverse magnetization out of plane My. 
A helical modulation with wavelength A ~ 85 /im is clearly 
visible. Simulation values: |c2|n3D — 320 Hz, Cduso ~ 0.8 Hz 
and q = 100 Hz (edges removed). 



with the wavelength of the helical modulation, for 
a derivation see Appendix C. In *^Rb with experimen- 
tally accessibly densities the wavelength is approximately 
SOfim and should be observable. Note that the scales for 
the two textures are related by S ~ A^(L/^s). 

Since the modulations of Fy and Fj_ decrease the total 
energy of those states, their regions of the phase diagram. 
Fig. [21 will be larger than predicted in the previous sub- 
section. However, the dipole strength must be large to 
introduce domains into the Fy state; and the energy gain 
in a helical texture relative to a uniform F± is small; 
so the phase boundaries will not change significantly at 
weak or moderate dipole strengths when we take these 
textures into account. 



V. NUMERICAL RESULTS 

We investigate numerically the ground state phase di- 
agram of a spin-1 condensate in external fields that give 
rise to a quadratic Zeeman shift and Larmor precession. 
The Metropolis algorithm^^ allows us to efficiently lo- 
cate minima of a given energy functional. We discretize 
the system on a lattice, and for the fundamental move we 
draw random deviations in the six real components of the 
field 5* from a normal distribution at a lattice site. The 
initial state is similarly generated from random normally 
distributed variables. 

A wide variety of simulation parameters {N, a, ay, 
Tmc, T^if^ M> Co, see below), for example 1 x 1 < 
A^ < 50 X 50, have been used to investigate the phase 
diagram(c2, q, Cd)- Energies have been calculated in Hz 
and the lengths have been inserted in /im. Unless other- 
wise noted, numerical results presented here use lattice 
constant a = 4 /im, thickness ay — 2 fim, and a system 
size of A^ = 30 X 30 plaquettes. We also add a chemi- 
cal potential to the energy, /i — 1202Hz/ini^^, in order 
to reproduce the experimental density for cq = l.9kHz. 
Finally, we set Tmc = 23 uK in the Metropolis weight 
g-{H) /kTMc ^ which strikes a good balance between re- 
ducing fluctuations and achieving convergence in a rea- 
sonable computation time and use a critical mean field 
temperature T^jp — IOOTmc- 

The phase diagram we have mapped out numerically 
agrees well with the results presented so far. In particu- 
lar, we have confirmed that the ferromagnetic states de- 
velop modulations governed by the strength of the dipole 
interaction. 

The algorithm described above tends to get trapped 
in local energy minima with varying densities of domain 
walls in the Fy region of the phase diagram. We can, 
however, locate the global minimum fairly confidently 
by starting the system in a variety of modulated states 
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(striped or checkerboard) and comparing the final en- 
ergies. The existence of metastable states as a conse- 
quence of dipolar interactions has been discussed before 
for spinor condensates in an optical lattice^^. We have 
not observed any tendencies for the simulation in the 
-Fj_/P|| region of the phase diagram to be trapped in a 
local energy minima, regardless of the initial configura- 
tion. This is as expected, since any possible local ground 
state configuration (Eq. ([7])) can smoothly turn into an- 
other, unlike in the case (Eq. ©). This symmetry 
between the two transverse components of the magneti- 
zation is present in the Hamiltonian without the dipole 
interaction, removed by the dipole interaction, and fi- 
nally restored by the rapid Larmor precession. However, 
even if the relaxational dynamics of the Metropolis al- 
gorithm used here does not apparently get trapped in 
a local minimum in this phase, the actual dynamics of 
the experimental system is primarily precessional rather 
than relaxational, which could lead to metastable states. 



A. Domain walls in 

Near the transition qc2{c2^ Cd), Eq. (HU, magnetization 
vortices with unit spin winding develop all the way along 
all domain walls, see Fig. [S] The vortices are alternating 
elliptical and hyperbolic Mermin-Ho vortices, with ferro- 
magnetic cores j^i^ The density of vortices increases with 
increasing dipole interaction, i.e. more domain walls ap- 
pear and the longitudinal length of each vortex decreases. 
The transverse length of the vortices increases with in- 
creasing quadratic Zeeman strength up to the transition 
line, which can be seen in the Fourier transform of the 
magnetization 



-irk^ 



M,{r, s). 



(27) 



as a rise in My{k^°'^); see Fig. [5] on the i^" side of the 
transition. The transition at qc2 itself remains sharp, and 
no vortices are observed for q > qc2- At a given instant 
in time does the perpendicular magnetization in all vor- 
tices in a domain boundary point in a specific direction. 
The correlations between the direction of the transverse 
magnetization of vortices in different domain walls are 
however weaker. 



B. Boundaries and trapping potential 

Finite size effects and the details of the trapping 
potential seem to have little impact on our results. The 
only finite size effect observed with hard-wall boundaries 
is a decrease in magnetization at the z — ±L/2 bound- 
aries in the transition from Fy to P±, as shown in Fig. |6l 
The approximative location of this transition line from 
the analytical calculation, Eq. (|16p . is |c2c|fi3i5 = 3.4Hz. 




30 32 



FIG. 5: (Color online) Transition to -F±/P|| from (a) 
For q slightly smaller than qc2, large Mermin-Ho vortices ap- 
pear between the stripes (plaquette size a — 4 /im). M^ix, z) 
is plotted on the horizontal axis, Mx{x,z) on the vertical 
axis and My{x,z) « for the whole region shown at this 
instant, (b) Consequently, the maximum value of the Fourier 
transform of the magnetization out of plane Afi(/s™°^), see 
Eq. (|27|) . increase before the phase transition. Simulation 
variables: |c2|n3D = 450 Hz, CdfizD ~ 7.2 Hz and q — 35 Hz 
(a), g = 30 - 39 Hz (b). 
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FIG. 6: (Color online) Transition to P± from F\\ . The par- 
allel magnetization Mz{z/a) is plotted for different values of 
|c2|?i3o = 2.4—4.8 Hz as a function of z/a. The magnetization 
is lowered at the boundaries around the transition point for 
a finite system. Simulation values: A'^ = 20 x 20, CdfL^D ~ 5.7 
Hz and q — —4 Hz. 
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We have also carried out simulations with an elliptical 
trap potential of the form [/(x) = U{vzi^)'^ + v^^^)'^), 
typically with U = 625 Hz/im~^ and Vz,Vx = 1 — 10 to 
more closely model experimental conditions. ^^"^ These 
simulations have shown no effect other than a decrease in 
the density and thereby related effects, as in the original 
paper of Ho on spinor condensates in optical traps. ^ 
For example, the wavelength of the helical modulation 
in F±/P\\ is inversely proportional to the density, see 
Fig. [7] which shows a change in wavelength through the 
condensate as the density changes. In particular, we 
have not seen the effect reported by Vengalattore et al^ 
in which the modulation wave vector is not aligned with 
the applied magnetic field but is instead influenced by 
the orientation of the trap. 



18 - 



5 10 15 20 25 30 

z/a 

FIG. 7: (Color online) Simulation of a helical modulated mag- 
netized condensate in an elliptical trap. n3D{x, z)Mz{x, z) is 
plotted on the horizontal axis, nsoix, z)Mx{x, z) on the ver- 
tical axis, and My(x,z) is a quarter of a wavelength ahead 
of Mx{x,z) as in Fig. ID but is not shown. The wavelength 
A(z) of the helical modulation increases with decreasing den- 
sity along the longitudinal axis. The distance between two 
neighbouring nodes is shown; the node to the left of them is 
outside the graph. Simulation parameters: = 1, Vx — 10, 
|c2|n3D = 540 Hz, Cdfi^D ~ 1-6 Hz and q — 120 Hz. 



VI. DISCUSSION 

We have mapped out the complete phase diagram for 
the model we have considered. Although the region oc- 
cupied by the phase moves and shrinks with the 
introduction of the dipole interaction, we find that it re- 
mains accessible at physical values of \c2\ and Cd in ^''Rb, 
for some values of the quadratic Zeeman shift q. Hence, 
by tuning q for ^''^Rb appropriately, the three phases F|| , 
F±/P\^ and P|| should be observable in experiments. We 
also find that a spatial modulations should be seen in at 
least the second of those phases. 

There are some disagreements between our result and 
other results obtained theoretically and more important 
experimentally. The length scale in the experiment is 
smaller than the pitch of the helical modulation we de- 
scribe above by a factor 10, roughly, for typical param- 
eters. Cherng and Demler— find a dynamical instability 
at a scale nearer that seen in experiment. That picture 



would suggest that even if the phase diagram obtained 
here describes the system at long times, the experimen- 
tal system might instead reach a long-lived metastable 
state. As explained in section V above, while we do see 
metastable states in some parts of the phase diagram, we 
do not see metastable checkerboard states in the region 
probed by current experiments, but this could be be- 
cause the Metropolis dynamics of our simulation is not 
the actual dynamics of the condensate, even if their ther- 
modynamics are the same. 

One challenge for this dynamical scenario is that in 
experiments, an imposed helical configuration with pitch 
A = 50 — 150/im2. quickly evolves into a state modulated 
at a smaller scale, again roughly ten times smaller than 
the stable, or at least metastable, supersolid state we pre- 
dict i^i^ This suggests that effects we have not taken into 
account prevent the current experimental system from 
finding this minimum. As an example, it is known that 
the dipole interaction makes the Larmor precession un- 
stable, according to Lamacraft^; as a result, the Larmor- 
averaged energy that is the main focus of the present 
work might not be an accurate description for long times. 

In order to observe the predicted supersolid clearly, 
our results suggest that the key is to suppress this Lar- 
mor instability while at the same time preserving the 
conservation of total magnetization in the field direction. 
The Larmor instability^ grows exponentially from ther- 
mal excitation of an initial perturbation at the Larmor 
frequency wl ■ Hence the time scale to reach a fixed final 
size of the instability is proportional to fvjjL/ [ksT) and 
can be increased either by increasing the magnetic field or 
decreasing the temperature. At the same time, an exper- 
iment should be designed to preserve the magnetization 
along the field direction for as long as possible, which re- 
quires a high degree of trap uniformity. One motivation 
for continued exploration of this system is that our re- 
sults show that the Larmor-averaged system does have a 
supersolid ground state for a wide range of parameters. 

Note added: As this work was being prepared for sub- 
mission, two e-prints appeared investigating the same 
experiment by slightly different approaches The 
first, by J. Zhang and T.-L. Ho, also investigates the 
static properties of ^^Rb using a deterministic numerical 
method and also gets the i^y state and a modulated F_\_ 
state. The main difference between their results and ours 
appears to be that they find a stripe phase rather than 
a helix for the phase with spins perpendicular to the ap- 
plied magnetic field. They find arrays of elliptical and hy- 
perbolic Mermin-Ho vortices, as a meta-stable dynamical 
state, between the stripes for the i^y state for all q. How- 
ever, they are smaller than the spin healing length and 
hence unobservable in our simulation, although we do 
see them close to the transition to the F^/ P\\ state. The 
second, by Y. Kawaguchi et al., finds a doubly periodic 
(checkerboard) spin pattern as a long-lived intermediate 
state through a combination of mean-field theory and 
numerical simulation of precession-averaged equations of 
motion. By adding energy dissipation to the dynamics. 
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they reach a stationary state similar to ours. 



Appendix A: The dipole term 

The dipolar energy of a magnetized fluid with magne- 
tization A4(x) is 



^ / dxdx' 



M ■ M' - 3{M ■ r){M' ■ r) 



where r = x — x' and A4' = Al(x'). The last term, or 
"s-wave" part, contributes to the contact interaction C2 
in the BEC Hamiltonian, and so should not be treated 
independently. In this paper we take the first, "d-wave" 
part to be the full dipolar interaction. This can, in turn, 
be decomposed into a "Coulomb" part that is positive 
semidefinite, and hence convenient for numerical work 
that searches for energy minima, and a contact part, as 
in Eq. (IT3l) . 

For both analytical and numerical work we need the di- 
mensionally reduced form of the Coulomb part expressed 
in a rotating frame. Ignoring the contact term in H^ip 
and performing two partial integrations we find 



Eg, .f/d3..3,. V-(.3.M(x))V.(.3.M(xO) ^ I / + 

I J d^xd^x'n,nMy{x,z)n2DMy{x',z') J dydy'MMM^ (Al) 

I 



where a{x,z) = dx{n2DMx{x, z)) + dz{n2DMz{x, z)) is 
an effective surface charge density. The density n2D has 
only a (x, z) dependence for a nonzero trapping potential 
U{x, z). The integrals over y can be performed explicitly 
for either Gaussian or boxcar profiles p; we choose the 
Gaussian form for the purposes of numerics. Then 

p{y)p{y') = -^e-^yl+y-^i< 
[dvP{y)][dy'P{y')] - ^^y+-f-\ -iyl+y'-y< (A2) 

with y± ^ y ±y' . The integrals over y^ are simple, and 
the integrals over y- can be put in terms of special func- 

^2 2 

tions with help of the identities / dx J^^^^.^ — e~ Ko{^) 

and / dx^^^=f — -^C/(i, 0, c^). Here Kq is a modified 
Bessel function and [/ is a confluent hypergeometric func- 
tion. 

For the numerics, discretize the remaining integrals as 
follows. Divide the 2-dimensional area into rectangular 
plaquettes and set the density n2D and magnetization M 
constant on each plaquette, 

M{x,z)^ M{a{r + ^),a{s + ^)), (A3) 

where a is the lattice constant and r, s are integers. Then 
do several variable substitutions. Going to variables x± 
and z± and scaling the coordinates by a allows us to 
replace 

/rp+l rq+1 
d'^xd'^x' -i' dx^ dz_(I-|a;_-p|)(I-|z_-(7|) 

(A4) 



since the integrands depend only on x_ , z_ . Here p — 
r' — r and q — s' — s. 

The integrals can then be computed numerically for 
< P,q < VN- The final step is to time-average the 
fields to take into account the rapid Larmor precession. 
This effectively means replacing 

a{p,q)a{p',q') ^ dz{n2DMz{p, q))dz>{n2DM^{p' ,q')) 
+ ]^dx{n2DMx{p,q))dx'{n2DMx{p' ,q')) 

+ ]^dx{n2DMy{p,q))dx'{n2DMy{p' ,q')) 

(A5) 

and 

My{p,q)My{p\q') ^ hl,{p,q)M,{p',q') 

1 (A6) 

+ -My{p,q)My{p',q') 

in Eq. (|AI[) . since the transverse components rotate into 
each other but the longitudinal component is unaffected. 



Appendix B: Helical modulation 

We can obtain a simple estimate of the wavelength 
of the transverse helical state to leading order in the 
strength of the dipole coupling by assuming a fully polar- 
ized time evolving state ^/'^(O, kzZ — jBot), Eq. with 
magnetization 

+ tMy - n2Dp(2/)e'("^^-''^«*). (BI) 
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Fourier transforming the kinetic and the dipole energy 
term, keeping only contributions that scale with the area 
of the two-dimensional system, the (areal) energy density 
of this state is 



energy 
area 



n2D 



k1 Cd n\jj f dky in 



2m 2 



2 2 



27r kl 



kl\~p{ky)? 
(B2) 

plus fc^ -independent terms. In the kinetic term, there is a 
factor of 1/2 because only half the atoms are in the mz = 
±1 states that carry kinetic energy. In the dipole term, 
the only extensive contribution to the energy comes from 
the out-of-plane component My, which gives a factor 1/2 
there as well. Notice also that the time dependence is 
gone. With k'^/{k^ + k: 
terms are 



1 - fc2/(P -f fcf), the relevant 



n2D 



2m 2 



22'^ 



du An 



27r 1 



■\~p{%\u)\\ (B3) 



and to lowest order in k^ we just need p(0) = jdy p{y) 
1 to arrive at 



n2D- — - TT— njclfel, 



■2m 2 "2 
which takes its minimum at 

TT n2DCd 



± 



2fi2/2 



m 



(B4) 



(B5) 



Appendix C: Dipole energy at uniform 
magnetization 

For a uniform condensate with maximal magnetiza- 
tion, aligned parallel to the magnetic field, the only con- 
tribution to the dipole energy comes from from the edges 
ai z — ±i/2. The second term in Eq. (|A1[) does not 
contribute and only the edges of the first 



In the limit L ^ W ^ T , the leading contribution to 
the energy comes solely from the first term, which de- 
scribes the self energy of two quasi-one-dimensional lines 
of charge. Indeed, it becomes just 



Egp = 2cdnla / dx^iW - x_)/x_ 

= 2cdnliyW\nW/T + 0(W) (C2) 

asymptotically, where the lower cutoff T has been chosen 
for convenience. 

The energy for the uniform out-of-plane configuration 

is 



E. 



dip 



(C3) 



Since there will be a term extensive in the planar size, it 
is simplest to ignore the effects of boundaries and work 
with a surface energy density 



2 o 



TT J dydy' J 



dr 



r[dypmdyp{y')] 



= 2ncdnlo / dydy'p{y)p{y')5{y - y') + 0(l/i?) 



— 2Trcdn2i 



T 



(C4) 



after integrating over the radial coordinate r followed by 
partial integration in y and y'. 
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